%{
AUTHOR: Felipe Arteaga
-------------------------------------------------------------------------
PROJECT: Warnings

El periodo de analisis es desde que se mando (
05-Sep-2020 23:35:58 UTC  para el resto) y que se manda el SMSImagen con mienduc
(07-Sep-2020 19:40:48, UTC). O sea una ventana de 1 dia y 20 horas
-------------------------------------------------------------------------
DESCRIPTION:
=========================================================================
%}


clearvars -except projectDir projectDirData fromMainWarningsPaper
clc;close all;fclose('all');feature('DefaultCharacterSet','UTF-8');

if(not(exist('projectDir','var')==1&&exist('projectDirData','var')==1&&exist('fromMainWarningsPaper','var')==1&&fromMainWarningsPaper))
    pcName=char(java.lang.System.getProperty('user.name'));
    if(strcmp(pcName,'felipe'))
        % PC Felipe
        myDir='/Users/felipe/Dropbox/';
        projectDir=[myDir,'git/warnings/'];
        projectDirData=[myDir,'projects/warnings/'];
        addpath(genpath([myDir,'/myMatlabFunctions/']));
    end
end
compileLatexTable=false;

%%
dirPlots=[projectDir,'/paper/figuresCL/RDs/'];
dirTable=[projectDir,'/paper/tablesCL/'];

saveFigFile=false;
withPlot=false;

%% Load data:
anho=2020;
anhoStr=sprintf('%i',anho);

% Load compiled admin/API/interventions data
load([projectDirData,'/data/chile/',anhoStr,'/inputRD'],'dataRD');


% if(ismember('DRiskExpostPreWEnd',dataRD.Properties.VariableNames))
%     error('Borrar esto');
% else
%     dataRD.DRiskExpostPreWEnd=dataRD.riskExpostEnd-dataRD.riskExpostPreW;
%     dataRD.DRiskExpostPreWPreSI=dataRD.riskExpostPreSI-dataRD.riskExpostPreW;
% end

% Make a shorter name, if not Stata complains
dataRD.treatedWarningWhatsapp_preSI=dataRD.treatedWarningWhatsapp_preSMSImage;
dataRD.placedInAnyPreferenceNWIE=dataRD.placedInAnyPreferenceNoWaitlistIniEnd;
dataRD.rural=dataRD.newMarket<0;

dataRD.totalTreatmentsPreSI=dataRD.totalPopups_IniPreSI+dataRD.treatedWarningWhatsapp+dataRD.treatedSMSPlain;
dataRD.totalTreatmentsEnd=dataRD.totalPopups+dataRD.treatedWarningWhatsapp+dataRD.treatedSMSPlain+dataRD.treatedSMSImage;

% Value added (this is from Chris)
vaAll=readtable([projectDirData,'/data/chile/auxiliar/FirmData.csv']);
va=vaAll(vaAll.Year==2016,:);% Last year on dataset
va.valueAddedEnrolled=va.VA2_AVE; % Not sure what is the difference between them
assert(allunique(va.School_RBD));
va.rbdEnrolled=va.School_RBD;
dataRD=outerjoin(dataRD,va,'keys',{'rbdEnrolled'},'mergeKeys',true,'type','left','rightVariables','valueAddedEnrolled');

data=dataRD;


%% RDs, OLSs and IVs
% Calls Stata to do all the inference, but plots and creates tables in
% matlab.

% Restricion to all specifications
% -Only applicants with some level of risk. Ommit mass points at extremes
% -Sample of the RCT (pobWhatsapp)
% -Avoid grade 9, experiment is not well implemented, many controls were
% treated
restrToAll=dataRD.riskWhatsapp>0.01&dataRD.riskWhatsapp<0.99&dataRD.pobWhatsapp==1&dataRD.grade<2;

subpobs=cell(1,3);
subpobs(1,:)={'Pooled',            restrToAll,'all'};
 %subpobs(2,:)={'w/cartilla',   restrToAll&(dataRD.intendedCartillaTWhatsapp==1|dataRD.controlWhatsapp==1),'wC'};
 %subpobs(3,:)={'wo/cartilla',      restrToAll&(dataRD.intendedCartillaCWhatsapp==1|dataRD.controlWhatsapp==1),'woC'};
 subpobs(2,:)={'w/image simple',   restrToAll&(dataRD.typeImageWhatsapp==1|dataRD.controlWhatsapp==1|dataRD.riskWhatsapp<=.3),'wI1'};
 subpobs(3,:)={'w/image bar',      restrToAll&(dataRD.typeImageWhatsapp==2|dataRD.controlWhatsapp==1|dataRD.riskWhatsapp<=.3),'wI2'};
 subpobs(4,:)={'w/image seats',    restrToAll&(dataRD.typeImageWhatsapp==3|dataRD.controlWhatsapp==1|dataRD.riskWhatsapp<=.3),'wI3'};
%subpobs(7,:)={'Control WhatsApp',restrToAll&dataRD.controlWhatsapp==1};



%% Define outcomes (and table order)
% Use 4rd row if want to start a new panel
% Use 5th row to define is IV is relevant. .5 if IV is not relevant, but
% plot is
% Use 6th row to define YLim of plot...optional!
% Use 7th row to define the enodgenous var in IV

outcomes={%'esSep','Economically Vulnerable','sep','A. Balance',0,[],'';...
    %'rural','Rural','rur','',0,[],'';...
    %'treatedWarningWhatsapp_preSI','WhatsApp read','fsPreSI','B. Message receipt',0,[],'';... % First stage of whatsapp intention over whatsapp read (moment right before SMSPlain was sent)
    %'treatedSMSImage','SMS reminder received','fs2','',0,[],'';... % First stage of whatsapp intention over SMSPlain.
    %'totalTreatmentsPreSI','Total treatments before final SMS','ttPreSP','',0.5,[-.04 6],'';...
    %'totalTreatmentsEnd','Total treatments endline','ttEnd','',0.5,[-.04 6],'';...
    %'changeAnyPreWPreSI','Any modification','amPreSMS','C. Outcomes in clean 44 hours before SMS followup',0.5,[-.04 .3],'';...
    'addAnyPreWPreSI','Add any (clean)','aaPreSMS','',0.5,[-.04 .2],'';...
    %'schoolsAddedPreWPreSI','Schools Added','saPreSMS','',1,[-.04 .3],'addAnyPreWPreSI';...
    %'riskExpostPreW','True Risk pre whatsapp','trPreW','',1,[-.04 .03]; ...
    %'riskExpostPreSI','True Risk post whatsapp','trPreSI','',1,[-.04 .03]; ...
    %'DRiskExpostPreWPreSI','$\Delta$ Risk','drPreSMS','',1,[-.04 .03],'addAnyPreWPreSI';...
    %'changeAnyPreWEnd','Any modification','amEnd','D. Endline outcomes',0.5,[-.04 .3],'';...
    'addAnyPreWEnd','Add any (endline)','aaEnd','',0.5,[-.04 .2],'';...
    %'schoolsAddedPreWEnd','Schools Added','saEnd','',1,[-.04 .3],'addAnyPreWEnd';...
    %'riskExpostEnd','True Risk final attempt','trEnd','',1,[-.04 .03]; ...
    %'DRiskExpostPreWEnd','$\Delta$ Risk','drEnd','',1,[-.04 .03],'addAnyPreWEnd'; ...
    %'assignedToPref','Placed to preference','ap','',1,[],'addAnyPreWEnd';...
    %'valueAddedEnrolled','Value Added Enrolled','vaE','',1,[],'addAnyPreWEnd';...
    %'placedInAddedIniEnd','Placed to added preference','apad','',1,[];...
    %'placedInAnyPreferenceNWIE','Placed to uncongested','apUncong','',1,[];...
    %'acceptOffer','Accept offer','ao','';...
    %'declineOffer','Decline offer','do','';...
    %'preferenciaAsign_1','Order assigned','oa','';...
    %'enrolledInAssigned','Enrolled in placed','ea','';...
    %'enrolledInAddedIniEnd','Enrolled in added','eaad',''...
    };


% Sample that has whatsapp (for RD on treated sample)
data.pobRD=data.intendedWhatsapp==1;

% Sample that has no whatsapp (for RD on control sample)
data.pobRD_onControl=data.controlWhatsapp==1;

% Sample that has warning (for RCT of whatsapp warning)
data.pobOLS=data.riskWhatsapp>.3;

% Sample that doesnt have warning (for RCT of wahtsapp plain, ie no warning)
data.pobOLS_notRisky=data.riskWhatsapp<.3;

% Restriction to all
inAnyPob=restrToAll;


commands=cell(size(subpobs,1)*size(outcomes,1),2);
counter=1;


% Special reg:
if(false)
    data.auxRest=restrToAll;
    stataExpress('rdrobust  DRiskExpostPreWPreSI riskWhatsapp if auxRest==1&pobRD==1, fuzzy(schoolsAddedPreWPreSI) c(.3) p(1) h(.1 .1)',data)
    stataExpress('rdrobust  DRiskExpostPreWEnd riskWhatsapp if auxRest==1&pobRD==1, fuzzy(schoolsAddedPreWEnd) c(.3) p(1) h(.1 .1)',data)
end


for p=1:size(subpobs,1)
    for r=1:2 % This loops over relevant or "placebo" populations (RCT on non risky, and RD on sample without whatsapp)
        if(r==1)
            % Relevants
            pOLS='pobOLS';
            pRD='pobRD';
            rL='';
            avoidIVs=false;
        elseif(r==2)
            % Placebo
            pOLS='pobOLS_notRisky';
            pRD='pobRD_onControl';
            rL='_placebo';
            avoidIVs=true;
        end
        
        
        data.(sprintf('subpop%i',p))=subpobs{p,2};
        
        %% OLSs and 2SLSs:
        
        for o=1:size(outcomes,1)
            
            % ITT
            commands(counter,:)={sprintf('reg_itt_%s_%i%s',outcomes{o,3},p,rL),...
                sprintf('reg  %s intendedWhatsapp if subpop%i==1&%s==1, robust',outcomes{o,1},p,pOLS)};counter=counter+1;
            
            if(outcomes{o,5}==1&&not(avoidIVs))
                % TOT
                
                commands(counter,:)={sprintf('reg_iv_%s_%i%s',outcomes{o,3},p,rL),...
                    sprintf('ivregress 2sls %s  (%s=intendedWhatsapp) if subpop%i==1&%s==1, first robust',outcomes{o,1},outcomes{o,7},p,pOLS)};counter=counter+1;
            end
        end
        
        %% RDs:
        
        
        for o=1:size(outcomes,1)
            % ITT
            commands(counter,:)={sprintf('rd_full_itt_%s_%i%s',outcomes{o,3},p,rL),...
                sprintf('rdrobust %s riskWhatsapp if subpop%i==1&%s==1, c(.3) p(2) h(.28 .68)',outcomes{o,1},p,pRD)};counter=counter+1;
            
            commands(counter,:)={sprintf('rd_local_itt_%s_%i%s',outcomes{o,3},p,rL),...
                sprintf('rdrobust %s riskWhatsapp if subpop%i==1&%s==1, c(.3) h(.1 .1)',outcomes{o,1},p,pRD)};counter=counter+1;
            
            
            if(outcomes{o,5}==1&&not(avoidIVs))
                % TOT
                commands(counter,:)={sprintf('rd_full_iv_%s_%i%s',outcomes{o,3},p,rL),...
                    sprintf('rdrobust %s riskWhatsapp if subpop%i==1&%s==1, c(.3) p(2) h(.28 .68) fuzzy(%s)',outcomes{o,1},p,pRD,outcomes{o,7})};counter=counter+1;
                commands(counter,:)={sprintf('rd_local_iv_%s_%i%s',outcomes{o,3},p,rL),...
                    sprintf('rdrobust %s riskWhatsapp if subpop%i==1&%s==1, c(.3) h(.1 .1) fuzzy(%s)',outcomes{o,1},p,pRD,outcomes{o,7})};counter=counter+1;
            end
        end
    end
end


% This is OK, since subsamples in each reagression depend on the label of
% the var, not in the position of the observations.
data=data(inAnyPob,:);

%% RUN all the regs and rds in Stata
res=stataCommand(commands,data,'selectcolumns',true);

%% Recover resutls and do tables/plots



cantPobs=(size(subpobs,1));
cantRegs=(size(outcomes,1));

% This is not "effective obs":
matrix_nl_t=nan(cantRegs,cantPobs);
matrix_nr_t=nan(cantRegs,cantPobs);

matrix_nl=nan(cantRegs,cantPobs);
matrix_nr=nan(cantRegs,cantPobs);
matrix_nT=nan(1,cantPobs);
matrix_nC=nan(1,cantPobs);

matrix_b=nan(cantRegs*2,cantPobs);
matrix_sd=nan(cantRegs*2,cantPobs);

matrix_pOrder=nan(cantRegs,cantPobs);
cell_pl=cell(cantRegs,cantPobs);
cell_pr=cell(cantRegs,cantPobs);


resStored=fieldnames(res);
close all


for p=1:cantPobs
         % N RCT
        matrix_nT(1,p)=sum(data.intendedWhatsapp&data.pobOLS&data.(sprintf('subpop%i',p)));
        matrix_nC(1,p)=sum((~data.intendedWhatsapp)&data.pobOLS&data.(sprintf('subpop%i',p)));
        
    for r=1:cantRegs
        
        withIv=false;
        
        
        res_reg_i=res.(sprintf('reg_itt_%s_%i',outcomes{r,3},p));
        res_rd_i=res.(sprintf('rd_local_itt_%s_%i',outcomes{r,3},p));
        
        if(p==1&&r==1)
            nRCT=mat2cellstr(res_reg_i.N,'rc',true);
        end
        
        % If is outcome (not first stage neither balance), then recover IV:
        if(outcomes{r,5}==1)
            withIv=true;
            resIv_reg_i=res.(sprintf('reg_iv_%s_%i',outcomes{r,3},p));
            resIv_rd_i=res.(sprintf('rd_local_iv_%s_%i',outcomes{r,3},p));
            
        end
        
        
        %% RCT
        % Beta
        matrix_b(1+(r-1)*2,p)=res_reg_i.b(1,1);
        % Sd Beta
        matrix_sd(1+(r-1)*2,p)=sqrt(res_reg_i.V(1,1));
   
        
        

        
        %% RD
        % Beta
        matrix_b(2+(r-1)*2,p)=res_rd_i.tau_cl;
        % Sd Beta
        matrix_sd(2+(r-1)*2,p)=res_rd_i.se_tau_cl;
        
        % N
        %         matrix_nl(r,p)=res_i.N_h_l;
        %         matrix_nr(r,p)=res_i.N_h_r;
        %
        %         matrix_nl_t(r,p)=res_i.N_l;
        %         matrix_nr_t(r,p)=res_i.N_r;
        
        % Bandwidth
        %matrix_nl(counter,p)=res.st(sirve);
        %matrix_nr(counter,p)=res.st(sirve);
        
        % Plot RD (+ control RD)
        if(withPlot)
            
            
            colorC=[0 0.4470 0.7410];
            colorT=[0.6350 0.0780 0.1840];
            % If I calculate IV means is an outcome, hence I wanna plot the
            % RD - RCT and not only the RD.
            plotRDControl=outcomes{r,5}>=.5;
            
            res_rd_i_full=res.(sprintf('rd_full_itt_%s_%i',outcomes{r,3},p));
            
            if(false&&plotRDControl)
                
                
                %%
                
                
                
                res_rd_i_control_full=res.(sprintf('rd_full_itt_%s_%i_placebo',outcomes{r,3},p));
                res_rd_i_control_local=res.(sprintf('rd_local_itt_%s_%i_placebo',outcomes{r,3},p));
                
                plotsubC=subpobs{p,2}&dataRD.controlWhatsapp==1;
                fC=plotRD(res_rd_i_control_full,dataRD,plotsubC,'color',colorC,'shiftPointEstimate',.15,'otherpointestimate',res_rd_i_control_local);
                ylim([-0.05 .15])
                xlabel('Predicted placement risk')
                ylabel(outcomes{r,2})
                easyExport([dirPlots,sprintf('%s_WhatsappCONTROL_%s_%s.png',anhoStr,outcomes{r,3},subpobs{p,3})]);
                close
                
                
                plotsubT=subpobs{p,2}&dataRD.controlWhatsapp==0;
                fT=plotRD(res_rd_i_full,dataRD,plotsubT,'color',colorT,'withPointEstimate',1,'otherpointestimate',res_rd_i);
                ylim([-0.05 .15])
                xlabel('Predicted placement risk')
                ylabel(outcomes{r,2})
                easyExport([dirPlots,sprintf('%s_WhatsappTREAT_%s_%s.png',anhoStr,outcomes{r,3},subpobs{p,3})]);
                close
                
                
            end
            
            
            
            
            figure
            if(plotRDControl)
                
                res_rd_i_control_full=res.(sprintf('rd_full_itt_%s_%i_placebo',outcomes{r,3},p));
                res_rd_i_control_local=res.(sprintf('rd_local_itt_%s_%i_placebo',outcomes{r,3},p));
                
                plotsubC=subpobs{p,2}&dataRD.controlWhatsapp==1;
                fC=plotRD(res_rd_i_control_full,dataRD,plotsubC,'color',colorC,'shiftPointEstimate',.15,'otherpointestimate',res_rd_i_control_local,'withcutoffverticalline',false,'marker','x','parameterLabel','\beta_{RD\_C}');
                hold on
            end
            
            plotsubT=subpobs{p,2}&dataRD.controlWhatsapp==0;
            fT=plotRD(res_rd_i_full,dataRD,plotsubT,'color',colorT,'withPointEstimate',1,'otherpointestimate',res_rd_i,'withcutoffverticalline',false,'parameterLabel','\beta_{RD\_T}');
            
            if(not(isempty(outcomes{r,6})))
                ylim(outcomes{r,6});
            end
            hold on;plot(.3*[1 1],ylim,'--','color',.5*[1 1 1]);hold off
            if(plotRDControl)
                
                % Add legend and ITT of RCT result:
                legend('Location','NorthEast');  % this sequence shows just the first 2 legends
                legend([fT.binscatter_l.scatter fC.binscatter_r.scatter],{'With WhatsApp (T)' 'Without WhatsApp (C)'},'Interpreter','latex');
                legend boxoff
                annotation('textbox',[.6 .65 .1 .1],'String',sprintf('$ITT_{RCT}: %.3f^{%s}$\n$\\quad\\quad\\quad\\quad(%.3f)$',res_reg_i.b(1,1),getStars(res_reg_i.b(1,1),sqrt(res_reg_i.V(1,1))),sqrt(res_reg_i.V(1,1))),'FitBoxToText','on','Interpreter','latex','edgecolor','none','color','k');
                
                hold off
                
            end
            
            if(not(isempty(outcomes{r,6})))
                ylim(outcomes{r,6});
            end
            
          
                
            
            
            xlabel('Predicted placement risk')
            ylabel(outcomes{r,2})
            easyExport([dirPlots,sprintf('%s_Whatsapp_%s_%s.png',anhoStr,outcomes{r,3},subpobs{p,3})]);
            
            if(saveFigFile)
                savefig(sprintf('%sW_%s_%s',dirPlots,outcomes{r,3},subpobs{p,3}))
            end
            
            
        end
    end
end





%%

posPan=find(not(cellfun(@isempty,outcomes(:,4))));
panel=[num2cell(posPan-1),outcomes(posPan,4)];



head={'RCT','RCT','RD','RD';'ITT','IV','ITT','IV'};
matTable=matrix_b;
matTable_sd=matrix_sd;
opt=struct;

%nRD=mat2cellstr(1000,'rc',true);
note=['ITT effects of 2020 WhatsApp warnings intervention. RCT columns: effects of random assignment to treatment group vs. control group for students with predicted risk $> 0.30$. Robust SEs in parentheses. N=',nRCT,'. '...
    'RD columns: regression discontinuity evaluation in treatment group around 0.30 cutoff. RD specifications computed using local linear fit with a bandwidth of 0.1. ',...
    'Standard errors are heteroskedasticity-robust nearest neighbor variance estimator with minimum of 3 neighbors, as in \citet{calonico2014robust}. ITT column shows effects of group assignment. IV columns show the instrumental variable specification, where the endogenous regressor is the add any school indicator, instrumented with group assignment for the RCT, and with a dummy of crossing the risky threshold for the RD. ',...
    'Panel A: balance tests on predetermined characteristics. Panel B: message receipt. ``WhatsApp read'''' is an indicator equal to one if applicant views the WhatsApp treatment message. ``SMS remainder received'''' is indicator for receiving SMS reminder 44 hours later. Panel C: outcomes within 44 hour window between WhatsApp intervention and followup SMS.  Panel D: endline choice behavior and placement outcomes. See section \ref{sec:RCT-eval} for details.  '];
opt.note=''; % note
opt.header=subpobs(:,1)';
assert(cantPobs==4)
opt.header={'Personalized','Personalized - By image type','Personalized - By image type','Personalized - By image type';'Pooled','Simple','Warning bar','Vacancy'};


%opt.stars=getStars(matTable,matTable_sd);

primeraCol=repmat({'','RCT';'','RD'},cantRegs,1);
primeraCol(1:2:end,1)=outcomes(:,2);

% Only RCT:
matTable=[matTable(1:2:end,:);matrix_nT;matrix_nC];
matTable_sd=[matTable_sd(1:2:end,:);nan(2,4)];
primeraCol=[primeraCol(1:2:end,1);'N Treatment';'N Control'];

opt.stderrs=mat2cellstr(matTable_sd,'conParentesis',true,'precisionDecimal','%.3f');
opt.primeracolumna=primeraCol;%[regressions(:,1);{'N'};regressions(:,1);{'NL';'NR'}];
opt.addColumnNumber=true;
opt.columnafantasma=1;
opt.filafantasma=cantRegs;
%opt.filafantasma=2:2:(cantRegs*2-1);
opt.adjust=true;
opt.title='RCT Estimates Split by WhatsApp Message';
opt.positionParameter='H';
opt.panel=panel;
opt.label='tab:2020RCTarms';
opt.file=sprintf('%s2020_tableArmsRCT',dirTable);

opt.note=' Evaluation of 2020 WhatsApp RCT, splitting out by message type. Treatment arms are as described in Online Appendix \ref{app:WhatsAppArms}. Online Appendix Figure \ref{fig:2020WhatsappImages} shows example images. Estimates are differences in the share of students adding any school to their baseline application between the treatment group and a control group that did not get any message.  Column 1 contains pooled estimates of the treatments from columns 2-4. ``Clean window'''' corresponds to outcomes measured 44 hrs after WhatsApp messages were sent. ``Endline'''' outcomes were measured at the end of the application process (75 hrs after WhatsApp message). See section \ref{sec:intdes} for a description of the WhatsApp RCT.';

tabla=cell2latex(mat2cellstr(matTable,'precisionDecimal','%.3f','revisarFilas',true),'opts',opt);
if(compileLatexTable)
compileLatex(tabla);
end